//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   testAngularMesh.cc
 * \author Jeremy Roberts
 * \date   06/18/2011
 * \brief  Test of AngularMesh class for several cases.
 * \note   Copyright (C) 2011 Jeremy Roberts
 */
//---------------------------------------------------------------------------//
// $Rev:: 120                                           $:Rev of last commit
// $Author:: j.alyn.roberts@gmail.com                   $:Author of last commit
// $Date:: 2011-07-14 13:55:43 +0000 (Thu, 14 Jul 2011) $:Date of last commit
//---------------------------------------------------------------------------//

#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cstdlib>


#include "angle/AngularMesh.hh"
#include "angle/Octant.hh"
#include "angle/QuadratureFactory.hh"
#include "angle/Definitions.hh"
#include "tran/Definitions.hh"
#include "tran/Traits.hh"
#include "../utilities/DBC.hh"
#include "../utilities/Soft_Equivalence.hh"


using namespace slabtran;
using namespace std;

int test_3d_level_symmetric(int order)
{
  //-----------------------------------------------------------------------//
  util::SP<AngularMesh<_3D> > d_octants;
  d_octants = new AngularMesh<_3D> (order, LEVELSYM, false);

  // array of verification values
  double v[24][4];
    v[ 0][0]= 0.3500211745815407; v[ 0][1]= 0.3500211745815407; v[ 0][2]= 0.8688903007222012; v[ 0][3]= 0.523598775598299;
    v[ 1][0]= 0.3500211745815407; v[ 1][1]= 0.8688903007222012; v[ 1][2]= 0.3500211745815406; v[ 1][3]= 0.523598775598299;
    v[ 2][0]= 0.8688903007222012; v[ 2][1]= 0.3500211745815407; v[ 2][2]= 0.3500211745815406; v[ 2][3]= 0.523598775598299;
    v[ 3][0]=-0.3500211745815407; v[ 3][1]= 0.3500211745815407; v[ 3][2]= 0.8688903007222012; v[ 3][3]= 0.523598775598299;
    v[ 4][0]=-0.3500211745815407; v[ 4][1]= 0.8688903007222012; v[ 4][2]= 0.3500211745815406; v[ 4][3]= 0.523598775598299;
    v[ 5][0]=-0.8688903007222012; v[ 5][1]= 0.3500211745815407; v[ 5][2]= 0.3500211745815406; v[ 5][3]= 0.523598775598299;
    v[ 6][0]= 0.3500211745815407; v[ 6][1]=-0.3500211745815407; v[ 6][2]= 0.8688903007222012; v[ 6][3]= 0.523598775598299;
    v[ 7][0]= 0.3500211745815407; v[ 7][1]=-0.8688903007222012; v[ 7][2]= 0.3500211745815406; v[ 7][3]= 0.523598775598299;
    v[ 8][0]= 0.8688903007222012; v[ 8][1]=-0.3500211745815407; v[ 8][2]= 0.3500211745815406; v[ 8][3]= 0.523598775598299;
    v[ 9][0]=-0.3500211745815407; v[ 9][1]=-0.3500211745815407; v[ 9][2]= 0.8688903007222012; v[ 9][3]= 0.523598775598299;
    v[10][0]=-0.3500211745815407; v[10][1]=-0.8688903007222012; v[10][2]= 0.3500211745815406; v[10][3]= 0.523598775598299;
    v[11][0]=-0.8688903007222012; v[11][1]=-0.3500211745815407; v[11][2]= 0.3500211745815406; v[11][3]= 0.523598775598299;
    v[12][0]= 0.3500211745815407; v[12][1]= 0.3500211745815407; v[12][2]=-0.8688903007222012; v[12][3]= 0.523598775598299;
    v[13][0]= 0.3500211745815407; v[13][1]= 0.8688903007222012; v[13][2]=-0.3500211745815406; v[13][3]= 0.523598775598299;
    v[14][0]= 0.8688903007222012; v[14][1]= 0.3500211745815407; v[14][2]=-0.3500211745815406; v[14][3]= 0.523598775598299;
    v[15][0]=-0.3500211745815407; v[15][1]= 0.3500211745815407; v[15][2]=-0.8688903007222012; v[15][3]= 0.523598775598299;
    v[16][0]=-0.3500211745815407; v[16][1]= 0.8688903007222012; v[16][2]=-0.3500211745815406; v[16][3]= 0.523598775598299;
    v[17][0]=-0.8688903007222012; v[17][1]= 0.3500211745815407; v[17][2]=-0.3500211745815406; v[17][3]= 0.523598775598299;
    v[18][0]= 0.3500211745815407; v[18][1]=-0.3500211745815407; v[18][2]=-0.8688903007222012; v[18][3]= 0.523598775598299;
    v[19][0]= 0.3500211745815407; v[19][1]=-0.8688903007222012; v[19][2]=-0.3500211745815406; v[19][3]= 0.523598775598299;
    v[20][0]= 0.8688903007222012; v[20][1]=-0.3500211745815407; v[20][2]=-0.3500211745815406; v[20][3]= 0.523598775598299;
    v[21][0]=-0.3500211745815407; v[21][1]=-0.3500211745815407; v[21][2]=-0.8688903007222012; v[21][3]= 0.523598775598299;
    v[22][0]=-0.3500211745815407; v[22][1]=-0.8688903007222012; v[22][2]=-0.3500211745815406; v[22][3]= 0.523598775598299;
    v[23][0]=-0.8688903007222012; v[23][1]=-0.3500211745815407; v[23][2]=-0.3500211745815406; v[23][3]= 0.523598775598299;
    // octant iterator
  AngularMesh<_3D>::iterator octant_itr = d_octants->begin();
  double sum_wt = 0.0;
  int idx = 0;
  for (; octant_itr != d_octants->end(); octant_itr++)
  {
    // get smart pointer to octant
    AngularMesh<_3D>::SP_octant octant = *octant_itr;
    // angle/weight iterators
    Octant<_3D>::angle_iterator angles = octant->angle_begin();
    Octant<_3D>::weight_iterator weights = octant->wt_begin();
    for (int m = 0; angles != octant->angle_end(); m++, angles++, weights++)
    {
      // compare using soft equivalence
      printf("%4i    %13.10f   %13.10f   %13.10f   %13.10f   \n", idx,
             (*angles)[MU], (*angles)[ETA], (*angles)[XI], *weights);
      if (!util::soft_equiv(v[idx][MU], (*angles)[MU])) return 1;
      if (!util::soft_equiv(v[idx][ETA], (*angles)[ETA])) return 2;
      if (!util::soft_equiv(v[idx][XI], (*angles)[XI])) return 3;
      if (!util::soft_equiv(v[idx][3], *weights)) return 4;
      sum_wt += *weights;
      ++idx;
    }
  }
  return 0;
}

int test_2d_level_symmetric(int order)
{
    //-----------------------------------------------------------------------//
    util::SP<AngularMesh<_2D> >    d_octants;
    d_octants = new AngularMesh<_2D>(order,LEVELSYM,false);

    // array of verification values
    double v[24][3];
    v[0][0] = 0.266635401516705;v[0][1] = 0.266635401516705;v[0][2] = 0.553316558825601;
    v[1][0] = 0.266635401516705;v[1][1] = 0.681507726536547;v[1][2] = 0.493880992370999;
    v[2][0] = 0.266635401516705;v[2][1] = 0.926180935517489;v[2][2] = 0.553316558825601;
    v[3][0] = 0.681507726536547;v[3][1] = 0.266635401516705;v[3][2] = 0.493880992370999;
    v[4][0] = 0.681507726536547;v[4][1] = 0.681507726536547;v[4][2] = 0.493880992370999;
    v[5][0] = 0.926180935517489;v[5][1] = 0.266635401516705;v[5][2] = 0.553316558825601;
    v[6][0] = -0.266635401516705;v[6][1] = 0.266635401516705;v[6][2] = 0.553316558825601;
    v[7][0] = -0.266635401516705;v[7][1] = 0.681507726536547;v[7][2] = 0.493880992370999;
    v[8][0] = -0.266635401516705;v[8][1] = 0.926180935517489;v[8][2] = 0.553316558825601;
    v[9][0] = -0.681507726536547;v[9][1] = 0.266635401516705;v[9][2] = 0.493880992370999;
    v[10][0] = -0.681507726536547;v[10][1] = 0.681507726536547;v[10][2] = 0.493880992370999;
    v[11][0] = -0.926180935517489;v[11][1] = 0.266635401516705;v[11][2] = 0.553316558825601;
    v[12][0] = 0.266635401516705;v[12][1] = -0.266635401516705;v[12][2] = 0.553316558825601;
    v[13][0] = 0.266635401516705;v[13][1] = -0.681507726536547;v[13][2] = 0.493880992370999;
    v[14][0] = 0.266635401516705;v[14][1] = -0.926180935517489;v[14][2] = 0.553316558825601;
    v[15][0] = 0.681507726536547;v[15][1] = -0.266635401516705;v[15][2] = 0.493880992370999;
    v[16][0] = 0.681507726536547;v[16][1] = -0.681507726536547;v[16][2] = 0.493880992370999;
    v[17][0] = 0.926180935517489;v[17][1] = -0.266635401516705;v[17][2] = 0.553316558825601;
    v[18][0] = -0.266635401516705;v[18][1] = -0.266635401516705;v[18][2] = 0.553316558825601;
    v[19][0] = -0.266635401516705;v[19][1] = -0.681507726536547;v[19][2] = 0.493880992370999;
    v[20][0] = -0.266635401516705;v[20][1] = -0.926180935517489;v[20][2] = 0.553316558825601;
    v[21][0] = -0.681507726536547;v[21][1] = -0.266635401516705;v[21][2] = 0.493880992370999;
    v[22][0] = -0.681507726536547;v[22][1] = -0.681507726536547;v[22][2] = 0.493880992370999;
    v[23][0] = -0.926180935517489;v[23][1] = -0.266635401516705;v[23][2] = 0.553316558825601;
    // octant iterator
    AngularMesh<_2D>::iterator octant_itr = d_octants->begin();
    double sum_wt = 0.0;
    int idx = 0;
    for (; octant_itr != d_octants->end(); octant_itr++)
    {
        // get smart pointer to octant
        AngularMesh<_2D>::SP_octant octant = *octant_itr;
        // angle/weight iterators
        Octant<_2D>::angle_iterator  angles  = octant->angle_begin();
        Octant<_2D>::weight_iterator weights = octant->wt_begin();
        for (int m = 0; angles != octant->angle_end(); m++, angles++, weights++)
        {
            //            printf ("%4i    %13.10f   %13.10f   %13.10f   \n",
            //                    idx, (*angles)[MU],(*angles)[ETA], *weights );
            // compare using soft equivalence
            if ( !util::soft_equiv( v[idx][MU], (*angles)[MU]  ) ) return 1;
            if ( !util::soft_equiv( v[idx][ETA],(*angles)[ETA] ) ) return 2;
            if ( !util::soft_equiv( v[idx][2], *weights        ) ) return 4;
            sum_wt += *weights;
            ++idx;
        }
    }
    return 0;
}

int test_3d_uniform_equal(int order)
{
    //-----------------------------------------------------------------------//
    util::SP<AngularMesh<_3D> >    d_octants;
    d_octants    = new AngularMesh<_3D>(order,UNIEQUAL,false);
    // array of verification values
    double v[32][4];
    v[0][0] = 0.894542510829527;v[0][1] = 0.370531640104871;v[0][2] = 0.250000000000000;v[0][3] = 0.392699081698724;
    v[1][0] = 0.370531640104871;v[1][1] = 0.894542510829527;v[1][2] = 0.250000000000000;v[1][3] = 0.392699081698724;
    v[2][0] = 0.611088871101869;v[2][1] = 0.253121298225658;v[2][2] = 0.750000000000000;v[2][3] = 0.392699081698724;
    v[3][0] = 0.253121298225659;v[3][1] = 0.611088871101869;v[3][2] = 0.750000000000000;v[3][3] = 0.392699081698724;
    v[4][0] = -0.894542510829527;v[4][1] = 0.370531640104871;v[4][2] = 0.250000000000000;v[4][3] = 0.392699081698724;
    v[5][0] = -0.370531640104871;v[5][1] = 0.894542510829527;v[5][2] = 0.250000000000000;v[5][3] = 0.392699081698724;
    v[6][0] = -0.611088871101869;v[6][1] = 0.253121298225658;v[6][2] = 0.750000000000000;v[6][3] = 0.392699081698724;
    v[7][0] = -0.253121298225659;v[7][1] = 0.611088871101869;v[7][2] = 0.750000000000000;v[7][3] = 0.392699081698724;
    v[8][0] = 0.894542510829527;v[8][1] = -0.370531640104871;v[8][2] = 0.250000000000000;v[8][3] = 0.392699081698724;
    v[9][0] = 0.370531640104871;v[9][1] = -0.894542510829527;v[9][2] = 0.250000000000000;v[9][3] = 0.392699081698724;
    v[10][0] = 0.611088871101869;v[10][1] = -0.253121298225658;v[10][2] = 0.750000000000000;v[10][3] = 0.392699081698724;
    v[11][0] = 0.253121298225659;v[11][1] = -0.611088871101869;v[11][2] = 0.750000000000000;v[11][3] = 0.392699081698724;
    v[12][0] = -0.894542510829527;v[12][1] = -0.370531640104871;v[12][2] = 0.250000000000000;v[12][3] = 0.392699081698724;
    v[13][0] = -0.370531640104871;v[13][1] = -0.894542510829527;v[13][2] = 0.250000000000000;v[13][3] = 0.392699081698724;
    v[14][0] = -0.611088871101869;v[14][1] = -0.253121298225658;v[14][2] = 0.750000000000000;v[14][3] = 0.392699081698724;
    v[15][0] = -0.253121298225659;v[15][1] = -0.611088871101869;v[15][2] = 0.750000000000000;v[15][3] = 0.392699081698724;
    v[16][0] = 0.894542510829527;v[16][1] = 0.370531640104871;v[16][2] = -0.250000000000000;v[16][3] = 0.392699081698724;
    v[17][0] = 0.370531640104871;v[17][1] = 0.894542510829527;v[17][2] = -0.250000000000000;v[17][3] = 0.392699081698724;
    v[18][0] = 0.611088871101869;v[18][1] = 0.253121298225658;v[18][2] = -0.750000000000000;v[18][3] = 0.392699081698724;
    v[19][0] = 0.253121298225659;v[19][1] = 0.611088871101869;v[19][2] = -0.750000000000000;v[19][3] = 0.392699081698724;
    v[20][0] = -0.894542510829527;v[20][1] = 0.370531640104871;v[20][2] = -0.250000000000000;v[20][3] = 0.392699081698724;
    v[21][0] = -0.370531640104871;v[21][1] = 0.894542510829527;v[21][2] = -0.250000000000000;v[21][3] = 0.392699081698724;
    v[22][0] = -0.611088871101869;v[22][1] = 0.253121298225658;v[22][2] = -0.750000000000000;v[22][3] = 0.392699081698724;
    v[23][0] = -0.253121298225659;v[23][1] = 0.611088871101869;v[23][2] = -0.750000000000000;v[23][3] = 0.392699081698724;
    v[24][0] = 0.894542510829527;v[24][1] = -0.370531640104871;v[24][2] = -0.250000000000000;v[24][3] = 0.392699081698724;
    v[25][0] = 0.370531640104871;v[25][1] = -0.894542510829527;v[25][2] = -0.250000000000000;v[25][3] = 0.392699081698724;
    v[26][0] = 0.611088871101869;v[26][1] = -0.253121298225658;v[26][2] = -0.750000000000000;v[26][3] = 0.392699081698724;
    v[27][0] = 0.253121298225659;v[27][1] = -0.611088871101869;v[27][2] = -0.750000000000000;v[27][3] = 0.392699081698724;
    v[28][0] = -0.894542510829527;v[28][1] = -0.370531640104871;v[28][2] = -0.250000000000000;v[28][3] = 0.392699081698724;
    v[29][0] = -0.370531640104871;v[29][1] = -0.894542510829527;v[29][2] = -0.250000000000000;v[29][3] = 0.392699081698724;
    v[30][0] = -0.611088871101869;v[30][1] = -0.253121298225658;v[30][2] = -0.750000000000000;v[30][3] = 0.392699081698724;
    v[31][0] = -0.253121298225659;v[31][1] = -0.611088871101869;v[31][2] = -0.750000000000000;v[31][3] = 0.392699081698724;

    // octant iterator
    AngularMesh<_3D>::iterator octant_itr = d_octants->begin();
    double sum_wt = 0.0;
    int idx = 0;
    for (; octant_itr != d_octants->end(); octant_itr++)
    {
        // get smart pointer to octant
        AngularMesh<_3D>::SP_octant octant = *octant_itr;
        // angle/weight iterators
        Octant<_3D>::angle_iterator  angles  = octant->angle_begin();
        Octant<_3D>::weight_iterator weights = octant->wt_begin();
        for (int m = 0; angles != octant->angle_end(); m++, angles++, weights++)
        {
            // compare using soft equivalence
            if ( !util::soft_equiv( v[idx][MU], (*angles)[MU]  ) ) return 1;
            if ( !util::soft_equiv( v[idx][ETA],(*angles)[ETA] ) ) return 2;
            if ( !util::soft_equiv( v[idx][XI], (*angles)[XI]  ) ) return 3;
            if ( !util::soft_equiv( v[idx][3], *weights        ) ) return 4;
            sum_wt += *weights;
            ++idx;
        }
    }
    return 0;
}

int test_2d_uniform_equal(int order)
{
    //-----------------------------------------------------------------------//
    util::SP<AngularMesh<_2D> > d_octants;
    d_octants = new AngularMesh<_2D>(order,UNIEQUAL,false);
    // array of verification values
    double v[36][3];
    v[0][0] = 0.952415708813758;v[0][1] = 0.255199020035362;v[0][2] = 0.349065850398866;
    v[1][0] = 0.697216688778397;v[1][1] = 0.697216688778396;v[1][2] = 0.349065850398866;
    v[2][0] = 0.255199020035362;v[2][1] = 0.952415708813758;v[2][2] = 0.349065850398866;
    v[3][0] = 0.836516303737808;v[3][1] = 0.224143868042013;v[3][2] = 0.349065850398866;
    v[4][0] = 0.612372435695795;v[4][1] = 0.612372435695794;v[4][2] = 0.349065850398866;
    v[5][0] = 0.224143868042014;v[5][1] = 0.836516303737808;v[5][2] = 0.349065850398866;
    v[6][0] = 0.533935590185808;v[6][1] = 0.143067610200522;v[6][2] = 0.349065850398866;
    v[7][0] = 0.390867979985286;v[7][1] = 0.390867979985286;v[7][2] = 0.349065850398866;
    v[8][0] = 0.143067610200522;v[8][1] = 0.533935590185808;v[8][2] = 0.349065850398866;
    v[9][0] = -0.952415708813758;v[9][1] = 0.255199020035362;v[9][2] = 0.349065850398866;
    v[10][0] = -0.697216688778397;v[10][1] = 0.697216688778396;v[10][2] = 0.349065850398866;
    v[11][0] = -0.255199020035362;v[11][1] = 0.952415708813758;v[11][2] = 0.349065850398866;
    v[12][0] = -0.836516303737808;v[12][1] = 0.224143868042013;v[12][2] = 0.349065850398866;
    v[13][0] = -0.612372435695795;v[13][1] = 0.612372435695794;v[13][2] = 0.349065850398866;
    v[14][0] = -0.224143868042014;v[14][1] = 0.836516303737808;v[14][2] = 0.349065850398866;
    v[15][0] = -0.533935590185808;v[15][1] = 0.143067610200522;v[15][2] = 0.349065850398866;
    v[16][0] = -0.390867979985286;v[16][1] = 0.390867979985286;v[16][2] = 0.349065850398866;
    v[17][0] = -0.143067610200522;v[17][1] = 0.533935590185808;v[17][2] = 0.349065850398866;
    v[18][0] = 0.952415708813758;v[18][1] = -0.255199020035362;v[18][2] = 0.349065850398866;
    v[19][0] = 0.697216688778397;v[19][1] = -0.697216688778396;v[19][2] = 0.349065850398866;
    v[20][0] = 0.255199020035362;v[20][1] = -0.952415708813758;v[20][2] = 0.349065850398866;
    v[21][0] = 0.836516303737808;v[21][1] = -0.224143868042013;v[21][2] = 0.349065850398866;
    v[22][0] = 0.612372435695795;v[22][1] = -0.612372435695794;v[22][2] = 0.349065850398866;
    v[23][0] = 0.224143868042014;v[23][1] = -0.836516303737808;v[23][2] = 0.349065850398866;
    v[24][0] = 0.533935590185808;v[24][1] = -0.143067610200522;v[24][2] = 0.349065850398866;
    v[25][0] = 0.390867979985286;v[25][1] = -0.390867979985286;v[25][2] = 0.349065850398866;
    v[26][0] = 0.143067610200522;v[26][1] = -0.533935590185808;v[26][2] = 0.349065850398866;
    v[27][0] = -0.952415708813758;v[27][1] = -0.255199020035362;v[27][2] = 0.349065850398866;
    v[28][0] = -0.697216688778397;v[28][1] = -0.697216688778396;v[28][2] = 0.349065850398866;
    v[29][0] = -0.255199020035362;v[29][1] = -0.952415708813758;v[29][2] = 0.349065850398866;
    v[30][0] = -0.836516303737808;v[30][1] = -0.224143868042013;v[30][2] = 0.349065850398866;
    v[31][0] = -0.612372435695795;v[31][1] = -0.612372435695794;v[31][2] = 0.349065850398866;
    v[32][0] = -0.224143868042014;v[32][1] = -0.836516303737808;v[32][2] = 0.349065850398866;
    v[33][0] = -0.533935590185808;v[33][1] = -0.143067610200522;v[33][2] = 0.349065850398866;
    v[34][0] = -0.390867979985286;v[34][1] = -0.390867979985286;v[34][2] = 0.349065850398866;
    v[35][0] = -0.143067610200522;v[35][1] = -0.533935590185808;v[35][2] = 0.349065850398866;
    // octant iterator
    AngularMesh<_2D>::iterator octant_itr = d_octants->begin();
    double sum_wt = 0.0;
    int idx = 0;
    for (; octant_itr != d_octants->end(); octant_itr++)
    {
        // get smart pointer to octant
        AngularMesh<_2D>::SP_octant octant = *octant_itr;
        // angle/weight iterators
        Octant<_2D>::angle_iterator  angles  = octant->angle_begin();
        Octant<_2D>::weight_iterator weights = octant->wt_begin();
        for (int m = 0; angles != octant->angle_end(); m++, angles++, weights++)
        {
            //            printf ("%4i    %13.10f   %13.10f   %13.10f   \n",
            //                    idx, (*angles)[MU],(*angles)[ETA], *weights );
            // compare using soft equivalence
            if ( !util::soft_equiv( v[idx][MU], (*angles)[MU]  ) ) return 1;
            if ( !util::soft_equiv( v[idx][ETA],(*angles)[ETA] ) ) return 2;
            if ( !util::soft_equiv( v[idx][2], *weights        ) ) return 4;
            sum_wt += *weights;
            ++idx;
        }
    }
    return 0;
}

int test_2d_quadruple(int order)
{
    //-----------------------------------------------------------------------//
    util::SP<AngularMesh<_2D> > d_octants;
    d_octants = new AngularMesh<_2D>(order,QUADRUPLE,false);

    // array of verification values
    double v[32][3];
    int i = 0;
    v[0+i][0]=4.46976567552685291442200e-02; v[0+i][1]=3.96864413885541210035759e-01; v[0+i][2]=1.39829431524808231013885e-01;
    v[1+i][0]=1.02344482002639025620492e-01; v[1+i][1]=9.08702733272681703905853e-01; v[1+i][2]=4.18269587013405891181035e-01;
    v[2+i][0]=1.99290928999251676505196e-01; v[2+i][1]=3.46095895882544646430440e-01; v[2+i][2]=2.53727177824340569108921e-01;
    v[3+i][0]=4.56317587472837769700214e-01; v[3+i][1]=7.92457765320398675821495e-01; v[3+i][2]=7.58970130432341927927481e-01;
    v[4+i][0]=3.46095895882544646430440e-01; v[4+i][1]=1.99290928999251676505196e-01; v[4+i][2]=2.53727177824340569108921e-01;
    v[5+i][0]=7.92457765320398675821495e-01; v[5+i][1]=4.56317587472837769700214e-01; v[5+i][2]=7.58970130432341927927481e-01;
    v[6+i][0]=3.96864413885541210035759e-01; v[6+i][1]=4.46976567552685291442200e-02; v[6+i][2]=1.39829431524808231013885e-01;
    v[7+i][0]=9.08702733272681703905853e-01; v[7+i][1]=1.02344482002639025620492e-01; v[7+i][2]=4.18269587013405891181035e-01;
    i = 8;
    v[0+i][0]=-4.46976567552685291442200e-02; v[0+i][1]=3.96864413885541210035759e-01; v[0+i][2]=1.39829431524808231013885e-01;
    v[1+i][0]=-1.02344482002639025620492e-01; v[1+i][1]=9.08702733272681703905853e-01; v[1+i][2]=4.18269587013405891181035e-01;
    v[2+i][0]=-1.99290928999251676505196e-01; v[2+i][1]=3.46095895882544646430440e-01; v[2+i][2]=2.53727177824340569108921e-01;
    v[3+i][0]=-4.56317587472837769700214e-01; v[3+i][1]=7.92457765320398675821495e-01; v[3+i][2]=7.58970130432341927927481e-01;
    v[4+i][0]=-3.46095895882544646430440e-01; v[4+i][1]=1.99290928999251676505196e-01; v[4+i][2]=2.53727177824340569108921e-01;
    v[5+i][0]=-7.92457765320398675821495e-01; v[5+i][1]=4.56317587472837769700214e-01; v[5+i][2]=7.58970130432341927927481e-01;
    v[6+i][0]=-3.96864413885541210035759e-01; v[6+i][1]=4.46976567552685291442200e-02; v[6+i][2]=1.39829431524808231013885e-01;
    v[7+i][0]=-9.08702733272681703905853e-01; v[7+i][1]=1.02344482002639025620492e-01; v[7+i][2]=4.18269587013405891181035e-01;
    i = 16;
    v[0+i][0]=4.46976567552685291442200e-02; v[0+i][1]=-3.96864413885541210035759e-01; v[0+i][2]=1.39829431524808231013885e-01;
    v[1+i][0]=1.02344482002639025620492e-01; v[1+i][1]=-9.08702733272681703905853e-01; v[1+i][2]=4.18269587013405891181035e-01;
    v[2+i][0]=1.99290928999251676505196e-01; v[2+i][1]=-3.46095895882544646430440e-01; v[2+i][2]=2.53727177824340569108921e-01;
    v[3+i][0]=4.56317587472837769700214e-01; v[3+i][1]=-7.92457765320398675821495e-01; v[3+i][2]=7.58970130432341927927481e-01;
    v[4+i][0]=3.46095895882544646430440e-01; v[4+i][1]=-1.99290928999251676505196e-01; v[4+i][2]=2.53727177824340569108921e-01;
    v[5+i][0]=7.92457765320398675821495e-01; v[5+i][1]=-4.56317587472837769700214e-01; v[5+i][2]=7.58970130432341927927481e-01;
    v[6+i][0]=3.96864413885541210035759e-01; v[6+i][1]=-4.46976567552685291442200e-02; v[6+i][2]=1.39829431524808231013885e-01;
    v[7+i][0]=9.08702733272681703905853e-01; v[7+i][1]=-1.02344482002639025620492e-01; v[7+i][2]=4.18269587013405891181035e-01;
    i = 24;
    v[0+i][0]=-4.46976567552685291442200e-02; v[0+i][1]=-3.96864413885541210035759e-01; v[0+i][2]=1.39829431524808231013885e-01;
    v[1+i][0]=-1.02344482002639025620492e-01; v[1+i][1]=-9.08702733272681703905853e-01; v[1+i][2]=4.18269587013405891181035e-01;
    v[2+i][0]=-1.99290928999251676505196e-01; v[2+i][1]=-3.46095895882544646430440e-01; v[2+i][2]=2.53727177824340569108921e-01;
    v[3+i][0]=-4.56317587472837769700214e-01; v[3+i][1]=-7.92457765320398675821495e-01; v[3+i][2]=7.58970130432341927927481e-01;
    v[4+i][0]=-3.46095895882544646430440e-01; v[4+i][1]=-1.99290928999251676505196e-01; v[4+i][2]=2.53727177824340569108921e-01;
    v[5+i][0]=-7.92457765320398675821495e-01; v[5+i][1]=-4.56317587472837769700214e-01; v[5+i][2]=7.58970130432341927927481e-01;
    v[6+i][0]=-3.96864413885541210035759e-01; v[6+i][1]=-4.46976567552685291442200e-02; v[6+i][2]=1.39829431524808231013885e-01;
    v[7+i][0]=-9.08702733272681703905853e-01; v[7+i][1]=-1.02344482002639025620492e-01; v[7+i][2]=4.18269587013405891181035e-01;
    // octant iterator
    AngularMesh<_2D>::iterator octant_itr = d_octants->begin();
    double sum_wt = 0.0;
    int idx = 0;
    for (; octant_itr != d_octants->end(); octant_itr++)
    {
        // get smart pointer to octant
        AngularMesh<_2D>::SP_octant octant = *octant_itr;
        // angle/weight iterators
        Octant<_2D>::angle_iterator  angles  = octant->angle_begin();
        Octant<_2D>::weight_iterator weights = octant->wt_begin();
        for (int m = 0; angles != octant->angle_end(); m++, angles++, weights++)
        {
            //            printf ("%4i    %13.10f   %13.10f   %13.10f   \n",
            //                    idx, (*angles)[MU],(*angles)[ETA], *weights );
            // compare using soft equivalence
            if ( !util::soft_equiv( v[idx][MU], (*angles)[MU]  ) ) return 1;
            if ( !util::soft_equiv( v[idx][ETA],(*angles)[ETA] ) ) return 2;
            if ( !util::soft_equiv( v[idx][2], *weights        ) ) return 4;
            sum_wt += *weights;
            ++idx;
        }
    }
    return 0;
}


int test_1d_gauss_legendre(int order)
{
  //-----------------------------------------------------------------------//
  util::SP<AngularMesh<_1D> > d_octants;
  d_octants = new AngularMesh<_1D> (order, GAUSSLEG, false);
  // array of verification values
  double v[8][2];
  v[0][0] = 0.9602898564975362316835609;
  v[1][0] = 0.7966664774136267395915539;
  v[2][0] = 0.5255324099163289858177390;
  v[3][0] = 0.1834346424956498049394761;
  v[4][0] = -0.9602898564975362316835609;
  v[5][0] = -0.7966664774136267395915539;
  v[6][0] = -0.5255324099163289858177390;
  v[7][0] = -0.1834346424956498049394761;
  v[0][1] = 0.1012285362903762591525314;
  v[1][1] = 0.2223810344533744705443560;
  v[2][1] = 0.3137066458778872873379622;
  v[3][1] = 0.3626837833783619829651504;
  v[7][1] = 0.3626837833783619829651504;
  v[6][1] = 0.3137066458778872873379622;
  v[5][1] = 0.2223810344533744705443560;
  v[4][1] = 0.1012285362903762591525314;
  // octant iterator
  AngularMesh<_1D>::iterator octant_itr = d_octants->begin();
  double sum_wt = 0.0;
  int idx = 0;
  for (; octant_itr != d_octants->end(); octant_itr++)
  {
    // get smart pointer to octant
    AngularMesh<_1D>::SP_octant octant = *octant_itr;
    // angle/weight iterators
    Octant<_1D>::angle_iterator angles = octant->angle_begin();
    Octant<_1D>::weight_iterator weights = octant->wt_begin();
    for (int m = 0; angles != octant->angle_end(); m++, angles++, weights++)
    {
      // compare using soft equivalence
      if (!util::soft_equiv(v[idx][MU], *angles)) return 1;
      if (!util::soft_equiv(v[idx][1], *weights)) return 4;
      sum_wt += *weights;
      ++idx;
    }
  }
  return 0;
}

// test driver; return nonzero for a failure; should work with CTest then.
int main(int argc, char *argv[])
{
  Require(argc==2);
  int test = atoi(argv[1]);
  cout << " ------ testAngularMesh ------ " << endl;

  int test_value = 1;
  int order = 1;

  if (test == 1 || test == 7)
  {
    //----------------------------------------------------------------
    // test 3-d level symmetric

    order = 4;
    cout << " TEST (3-d Level Symmetric, N = " << order << "): ";
    test_value = test_3d_level_symmetric(order);
    if (test_value == 0) cout << " PASSED." << endl;
    else
      cout << " FAILED." << test_value << endl;
  }
  if (test == 2 || test == 7)
  {
    //----------------------------------------------------------------
    // test 2-d level symmetric
    order = 6;
    cout << " TEST (2-d Level Symmetric, N = " << order << "): ";
    test_value = test_2d_level_symmetric(order);
    if (test_value == 0) cout << " PASSED." << endl;
    else
      cout << " FAILED. " << test_value << endl;
  }
  if (test == 3 || test == 7)
  {
    //----------------------------------------------------------------
    // test 3-d uniform equal weight
    order = 4;
    cout << " TEST (3-d Uniform Equal, N = " << order << "): ";
    test_value = test_3d_uniform_equal(order);
    if (test_value == 0) cout << " PASSED." << endl;
    else
      cout << " FAILED." << test_value << endl;
  }
  if (test == 4 || test == 7)
  {
    //----------------------------------------------------------------
    // test 2-d uniform equal weight
    order = 6;
    cout << " TEST (2-d Uniform Equal, N = " << order << "): ";
    test_value = test_2d_uniform_equal(order);
    if (test_value == 0) cout << " PASSED." << endl;
    else
      cout << " FAILED. " << test_value << endl;
  }
  if (test == 5 || test == 7)
  {
    //----------------------------------------------------------------
    // test 2-d symmetric quadruple-range
    order = 4;
    cout << " TEST (2-d Symmetric Quadruple Range, N = " << order << "): ";
    test_value = test_2d_quadruple(order);
    if (test_value == 0) cout << " PASSED." << endl;
    else
      cout << " FAILED. " << test_value << endl;
  }
  if (test == 6 || test == 7)
  {
    //----------------------------------------------------------------
    // test 1-d gauss legendre
    order = 8;
    cout << " TEST (1-d Gauss-Legendre, N = " << order << "): ";
    test_value = test_1d_gauss_legendre(order);
    if (test_value == 0) cout << " PASSED." << endl;
    else
      cout << " FAILED." << test_value << endl;
  }
  if (test > 7 || test < 1)
  {
    cout << " test " << test << " not available " << endl;
    return 1;
  }

  Check(8==8);

  return test_value;
}
